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Abstract This paper is concerned with approximations of the Boltzmann equation based on the 
method of moments. We propose a generalization of the setting of the moment-closure problem 
from relative entropy to (^-divergences and a corresponding closure procedure based on minimiza¬ 
tion of (^-divergences. The proposed description encapsulates as special cases Grad’s classical closure 
based on expansion in Hermite polynomials and Levermore’s entropy-based closure. We establish 
that the generalization to divergence-based closures enables the construction of extended thermo¬ 
dynamic theories that avoid essential limitations of the standard moment-closure formulations such 
as inadmissibility of the approximate phase-space distribution, potential loss of hyperbolicity and 
singularity of flux functions at local equilibrium. The divergence-based closure leads to a hierarchy 
of tractable symmetric hyperbolic systems that retain the fundamental structural properties of the 
Boltzmann equation. 

Keywords Boltzmann equation • kinetic theory • moment closure • hyperbolic systems • entropy • 
divergence 


1 Introduction 

The Boltzmann equation provides a description of the molecular dynamics of fluid flows based on their 
one-particle phase-space distribution. However, the Boltzmann equation also encapsulates all con¬ 
ventional macroscopic flow models in the sense that its limit solutions correspond to solutions of the 
compressible Euler and Navier-Stokes equations mm- the incompressible Euler and Navier-Stokes 
equations [181135] , the incompressible Stokes equations [36] and the incompressible Navier-Stokes- 
Fourier system [ST]; see m for an overview. The Boltzmann equation is uniquely suited to describe 
flows in the transitional molecular/continuum regime and the corresponding rarefaction effects, by 
virtue of its inherent characterization of deviations of the velocity distribution from local equilibrium. 
Applications in which rarefaction effects play a significant role are multitudinous, including gas flow 
problems involving large mean free paths in liigh-altitude flows and hypobaric applications such as 
chemical vapor deposition; see mm and references therein for further examples. Moreover, the 
perpetual trend toward miniaturization in science and technology renders accurate descriptions of 
fluid flows in the transitional molecular/continuum regime of fundamental technological relevance, 
for instance, in nanoscale applications, micro-channel flows or flow in porous media [21]. The Boltz¬ 
mann equation also provides a prototype for kinetic models in many other applications that require 
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a description of the collective behavior of large ensembles of small particles, for instance, in semi¬ 
conductors [25], in plasmas and fusion and fission devices m and in dispersed-particle flows such as 
in fluidized-bed reactors mmm- 

Numerical approximation of the Boltzmann equation poses a formidable challenge, on account of 
its high dimensional setting: for a problem in D spatial dimensions, the one-particle phase-space is 2D 
dimensional. The corresponding computational complexity of conventional discretization methods for 
(integro-)differential equations, such as finite-element methods with uniform meshes, is prohibitive. 
Numerical approximations of the Boltzmann equations have been predominantly based on particle 
methods, such as the Direct Simulation Monte Carlo (DSMC) method [5][6]. Convergence proofs 
for these methods [56] however convey that their computational complexity depends sensitively on 
the Knudsen number, and the computational cost becomes prohibitive in the fluid-dynamical limit. 
Moreover, from an approximation perspective, DSMC can be inefficient, because it is inherent to the 
underlying Monte-Carlo process that the approximation error decays only as n ~2 as the number of 
simulation molecules, n, increases; see, for instance, m Thm. 5.14], Efficient computational model¬ 
ing of fluid flows in the transitional molecular/continuum regime therefore remains an outstanding 
challenge. 

An alternative approximation technique for the Boltzmann equation, which has been only rela¬ 
tively sparsely investigated, is the method of moments mmm- The method of moments represents 
a general statistical approximation technique which identifies parameters of an approximate distri¬ 
bution based on its moments [38]. Application of the method of moments to the Boltzmann equation 
engenders an evolution equation for the moments (weighted averages) of the phase-space distribution. 
An approximation based on moments is generally consistent with a restricted interest in functionals 
of the distribution corresponding to macroscopic properties of the fluid. The method of moments is 
closely related to extended thermodynamics; see, for instance, mm- 

Intrinsic to the method of moments is an approximation of the moment-closure relation, viz., a 
relation that closes the evolution equation for the moments. Moment-closure approximations for the 
Boltzmann equation were originally conceived by Grad [T9]. Grad’s moment closure is based on an 
expansion of the one-particle distribution in Hermite polynomials. For a linear Boltzmann equation 
extended with an exogenous forcing, Schmeiser and Zwirchmayr [45] have shown that the distribu¬ 
tion in Grad’s moment equations converges to the distribution of the underlying kinetic model as 
the order of the moment approximation tends to infinity, and to the solution of a corresponding 
drift-diffusion model in the macroscopic limit, i.e. as the Knudsen number tends to zero. Grad’s 
moment systems are impaired by two essential deficiencies, however, viz., the potential occurrence 
of inadmissible locally negative phase-space distributions and potential loss of hyperbolicity nnusii- 
Levermore [33] has developed a moment-closure procedure based on constrained entropy minimiza¬ 
tion, similar to Dreyer’s maximum-entropy closure in extended thermodynamics [15] • The entropy 
minimization procedure formally leads to an exponential closure. Levermore’s moment systems re¬ 
tain the fundamental structural properties of the Boltzmann equation, viz., conservation of mass, 
momentum and energy, Galilean invariance and entropy dissipation. Moreover, the moment systems 
form a hierarchy of symmetric hyperbolic systems and the corresponding distributions are non¬ 
negative. It was later shown by Junk [26], however, that Levermore’s moment-closure procedure is 
impaired by a realizability problem, in that there exist moments for which the minimum-entropy 
distribution is non-existent. Moreover, the fluxes in Levermore’s moment systems may become ar¬ 
bitrarily large in the vicinity of (local) equilibrium. Recent results by Junk [27], Schneider [49] and 
Pavan [43] convey that potential non-existence of solutions to the entropy-minimization problem can 
be avoided by relaxing the constraints to allow inequalities in the liighest-order moments. Moreover, 
the solution to the relaxed entropy minimization problem coincides with the solution to the origi¬ 
nal constrained entropy minimization problem if the latter admits a solution. The relaxation of the 
entropy-minimization problem however generally engenders the loss of a one-to-one correspondence 
between the moments and the distribution. Moreover, relaxation of the entropy minimization prob¬ 
lem does not resolve the potential singularities in the flux function, as these singularities are intrinsic 
to the exponential form of the closure relation. Another fundamental complication, pertaining to the 
implementation of moment systems based on exponential closure, is that the resulting formulation 
requires the evaluation of moments of exponentials of polynomials of, in principle, arbitrary order. It 
is generally accepted that the derivation of closed-form expressions for such moments is intractable, 
and accurate approximation of the moments is a notoriously difficult problem; see, for instance, [32!. 
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In this paper we consider alternative moment-closure relations for the Boltzmann equation, based 
on approximations of the exponential function derived from truncations of its standard limit defini¬ 
tion, exp(-) = lim n _ ) . 00 (l + ■/n) n . It is to be noted that closure relations derived from a series- 
expansion definition of the exponential have received scant attention before, e.g., by Brini and 
Ruggeri HD. Our motivation for considering the limit definition instead of the series-expansion 
definition for constructing the moment closures is based on the direct availability of a corresponding 
inverse relation for higher order approximations. We propose a generalization of the setting of the 
moment-closure problem from Kullback-Leibler divergence m (i.e relative entropy) to the class 
of (^-divergences [13] • The considered (^-divergences constitute an approximation to the Kullback- 
Leibler divergence in the vicinity of some Maxwellian. It will be shown that the approximate- 
exponential closure relation can be derived via constrained minimization of a corresponding (/^-divergence. 
The proposed description encapsulates as special cases Grad’s closure relation and Levermore’s 
entropy-based closure relation. For even order approximations of the exponential, the closure relation 
engenders non-negative phase-space distributions. Moreover, the corresponding moment systems are 
symmetric hyperbolic and tractable, in the sense that the formulation only requires the evaluation 
of higher-order moments of Gaussian distributions. The moment systems furthermore dissipate an 
appropriate (^-divergence, analogous to the dissipation of relative entropy of the Boltzmann equation, 
provided that the collision operator dissipates the corresponding (^-divergence. We will show that the 
class of collision operators that dissipate appropriate (/3-divergences includes the standard BGK [3] 
and generalized BGK [33] operators. 

The remainder of this paper is organized as follows. Section [2] abstracts, for completeness, well 
known structural features of the Boltzmann equation to be retained in the developed moment system 
approximation. Section [3] introduces concepts relevant to moment systems pertaining to subspace 
approximations and reviews the moment closures of Grad m and Levermore [33] in light of the 
aforementioned issues, namely, admissibility of phase-space distributions, hyperbolicity, realizability 
and tractability. Section [4] presents a novel tractable moment closure approximation and, moreover, 
it will be shown that the corresponding closed system of moment equations are well-posed and retain 
the structural features of the Boltzmann equation. Finally, section [6] gives a concluding discussion. 


2 The Boltzmann Equation 

Consider a monatomic gas, i.e. a gas composed of a single species of identical classical particles, 
contained in a fixed spatial domain Q C . Kinetic theory describes the state of such a gas by a 
non-negative (phase-space) density / = f(t,x,v) over the single-particle phase space fi x R . The 
evolution of / is considered to be governed by the Boltzmann equation, 

dtf + v i d x J=C(f) (1) 

where the collision operator / C(f) acts only on the v — (v\,... ,vn) dependence of / locally 
at each (t, x ) and the summation convention applies to repeated indices. The collision operator is 
assumed to possess certain conservation, symmetry and dissipation properties, viz., conservation of 
mass, momentum and energy, invariance under Galilean transformations and dissipation of appropri¬ 
ate entropy functionals. These fundamental properties of the collision operator are treated in further 
detail below. Our treatment of the conservation and symmetry properties is standard (see, for in¬ 
stance, [3D) and is presented merely for coherence and completeness. For the entropy-dissipation 
property, we consider a generalization of the usual (relative) entropy to (^-divergences [13] , to enable 
an exploration of the moment-closure problem in an extended setting; see Section [IJ 

To elaborate the conservation properties of the collision operator, let (•) denote integration in the 
velocity dependence of any scalar, vector or matrix valued measurable function over D-dimensional 
Lebesgue measure. A function x/j : R D —» R is called a collision invariant of C if 

<W)>= 0 v/e^(C), (2) 

where 0(C) C L^R , R>o) denotes the domain of C, which we consider to be a subset of the 
almost everywhere nonnegative Lebesgue integrable functions on R 15 . Equation Jl]) associates a 
scalar conservation law with each collision invariant: 

+d Xi (vii>f) =0 


(3) 
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We insist that { 1 ,Ui,...,ud,|u| 2 } are collision invariants of C and that the span of this set contains 
all collision invariants, i.e. 

(ipC(f))=0 V/ G @{C) ^ ip G spanjl, Vi,. .., vp, |u| 2 } =: E. 

The moments (/), (vif) and ( \v\ 2 f ), correspond to mass-density, the (components of) momentum- 
density and energy-density, respectively. Accordingly, the conservation law © implies that dU) con¬ 
serves mass, momentum and energy. 

The assumed symmetry properties of the collision operator pertain to commutation with trans¬ 
lational and rotational transformations. In particular, for all vectors u G and all orthogonal 
tensors O : R 13 — » R , we define the translation transformation T u : S>(C) —> @{C) and the rotation 
transformation To '■ @{C) —> S>{C) by: 

(T u f)(v) = f(u-v) Vfe9(C) 

(Tof)(v) = f(0*v) Vfe®(C) 

with O* the Euclidean adjoint of O. Note that the above transformations act on the v-dependence 
only. It is assumed that C possesses the following symmetries: 


C(T u f) = T u C(f), C(Tof) = ToC(f) (4) 

The symmetries Q) imply that (jT]) complies with Galilean invariance, i.e. if f(t,x,v) satisfies the 
Boltzmann equation ©> then for arbitrary u G R 13 and arbitrary orthogonal O : M 13 —>• M. D , so do 
f(t, x — ut , v — u) and f(t, 0*x, 0*v). 

The entropy dissipation property of C is considered in the extended setting of |33| Sec. 7], from 
which we derive the following definition: a convex function 77 : M>o —> M is called an entropy density 
for C if 

W(f)C(f))< 0 , V/ e 3>(C) (5) 

with rf(f) the derivative of ??(/), and if for every / G S>{C) the following equivalences hold: 

C(f) — 0 ^ (rf(f)C(f))= 0 ^ V \f)e E ( 6 ) 

Relation © implies that C dissipates the local entropy which leads to an abstraction of 

Boltzmann’s H-theorem for ©, asserting that solutions of the Boltzmann equation © satisfy the 
local entropy-dissipation law: 

dt(v(f))+d Xi (v iV (f)) = (C(f)v'(f)) <0. ( 7 ) 


The functions (??(/)), (viV(f)) an d (?/(/)£(/)} are referred to as entropy density, entropy flux and 
entropy-dissipation rate, respectively. The first equivalence in ([ 6 ]) characterizes local equilibria of C by 
vanishing entropy dissipation, while the second equivalence indicates the form of such local equilibria. 
For spatially homogeneous initial data, /o, Equations © and ([ 6 ]) suggest that equilibrium solutions, 
/eq, of (JTJ) are determined by: 


fe q = arg min {(??(/)} : / G @(C),{fip) = </oV’)}, (8) 

Equation Q identifies equilibria as minimizersQ of the entropy, subject to the constraint that the 
invariant moments are identical to the invariant moments of the initial distribution. 

The standard definition of entropy corresponds to a density / 1 —> /log/, possibly augmented 
with ftp where ip G E is any collision invariant. It is to be noted that for Maxwellians A4, i.e. 
distributions of the form 


M{v) := M^ UtT) {v) 


Q 

(27 rRT)^ 


exp 


V — u \ 2 \ 

2 RT ) 


(9) 


for some (g,u,T) G R>o x 1R 13 x R>o and a certain gas constant R G R>o, it holds that logA4 G E. 
Therefore, the relative entropy (/log (//A!)} of / with respect to M. is equivalent to (/log/) in 
the sense of dissipation characteristics. The physical interpretation of the entropy (/log/), due to 


1 We adopt the sign convention of diminishing entropy. 
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Boltzmann [7118119]. is that of a measure of degeneracy of macroscopic states, i.e. of the number of 
microscopic states that are consistent with the macroscopic state as described by the one-particle 
marginal, /. In the context of information theory, Shannon [50] showed that for discrete probability 
distributions, the density / i-> / log / is uniquely defined by the postulates of continuity, strong ad¬ 
ditivity and the property that mr\(\jrn) < nr/(l/n) whenever n < m. These postulates ensure that 
for discrete probability distributions the entropy yields a meaningful characterization of information 
content and, accordingly, rationalize an interpretation of entropy as a measure of the uncertainty or, 
conversely, information gain pertaining to an observation represented by the corresponding proba¬ 
bility distribution [23]. Kullback and Leibler [31] generalized Shannon’s definition of information to 
the abstract case and identified the divergencqj 


-Dkl(mi|M2) 


/i l°g(/i// 2 ) dv 


( 10 ) 


as a distance between mutually absolutely continuous measures p,\ and both absolutely continuous 
with respect to the measure v with Radon-Nikodym derivatives /i = dpi/dv and / 2 = dp, 2 /dv. 
The Kullback-Leibler divergence characterizes the mean information for discrimination between pi 
and /i 2 per observation from Noting that the Kullback-Leibler divergence m coincides with the 
relative entropy of /i with respect to / 2 , the relative entropy (/ log(//.M)) can thus be understood 
as a particular measure of the divergence of the one-particle marginal relative to the reference (or 
background ) distribution JA. Kullback-Leibler divergence was further generalized by Csiszar [T3] and 
Ali et. al. [I], who introduced a general class of distances between probability measures, referred to 
as ip- divergences, of the form: 

D v (n 1 \n 2 ) = J hy{h/f2)dv (li) 

where ip is some convex function subject to <p{ 1) = tp'(l) = 0 and tp"( 1) > 0. Note that the Kullback- 
Leibler divergence corresponds to the specific case ^kl(-) = (•) log(-). 

In this work, we depart from the standard (relative) entropy for (JTJ) and instead consider entropies 
based on particular (/^-divergences. These (^-divergences generally preclude the usual physical and 
information-theoretical interpretations, but still provide a meaningful entropy density in accordance 
with ([5J and ([6]). The considered (^-divergences yield a setting in which entropy-minimization based 
moment-closure approximations to © are not impaired by non-realizability, exhibit bounded fluxes 
in the vicinity of equilibrium, and are numerically tractable. 

Remark 1 Implicit to our adoption of (/?-divergence-based entropies is the assumption that such 
entropies comply with © and dn]) for a meaningful class of collision operators. It can be shown that 
the class of admissible collision operators includes the BGK operator [4] : 

Cbgk(/) = — S f ) (12) 

where r E R>o is a relaxation time and £(.) corresponds to the map /o (->■ / eq defined by ([8]). 
The Kuhn-Tucker optimality conditions associated with © convey that f]'(£f) E E and, therefore, 
(?/(£/)(/ ~ £/)) = 0. The dissipation inequality © then follows from the convexity of r/(-): 

(v'(f)CBG k(/)} = -r _1 ((??'(/) -?/(£/))(/-£/)) < 0 (13) 

Moreover, because equality in m holds if and only if / = £f, the condition {rf{f) Cbgk(/)) = 0 
implies that f — £f, which in turn yields Cbgk(/) = 0 and r/'(f) E E. The equivalences in © are 
therefore also verified. A similar result holds for the multi-scale generalization of the BGK operator 
introduced in [33]; see Appendix [Al 


2 The conventional definition of Kullback-Leibler divergence according to is historically incorrect, as 

Kullback and Leibler in fact referred to the symmetrization of (nsj as the “divergence”. 
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3 Moment Systems 

Moment systems are approximations of the Boltzmann equation based on a finite number of velocity- 
moments of the one-particle marginal. An inherit aspect of moment equations derived from m is 
that low-order moments are generally coupled with higher-order ones, and consequently a closed set 
of equations for the moments cannot be readily formulated. Therefore, a closure relation is required. 

To derive the moment equations from m and elaborate on the corresponding moment-closure 
problem, let M denote a finite-dimensional subspace of D- variate polynomials and let {mi(v)}^ =1 
represent a corresponding basis. Denoting the column M-vector of these basis elements by m, it 
holds that the moments {(mif)}f£ 1 of the one-particle marginal satisfy: 

d t (mf) +d Xi (vimf) = ( mC(f)) (14) 

It is to be noted that we implicitly assume in m that / resides in 

F := {/ € 3>(C) : mf G L 1 (R d ), vmf G L 1 (M £) ,M 15 ), mC{f) G T 1 (M' D ) for all m G M} (15) 

almost everywhere in the considered time interval (0, T) and the spatial domain 12. This assumption 
has been confirmed in specific settings of m but not for the general case; see [33] Sec. 4] and 
the references therein for further details. The moment-closure problem pertains to the fact that m 
provides only M relations between (2 + D)M independent variables, viz., the densities ( rtiif ), the flux 
components (viiriif) and the production terms ( rriiC(f )). Therefore, (1 + D)M auxiliary relations 
must be specified to close the system. Generally, moment systems are closed by expressing the 
fluxes and production terms as a function of the densities. Moment systems are generally closed by 
constructing an approximation to the distribution function from the densities and then evaluating 
the fluxes and production terms for the approximate distribution. Denoting by A C K M a suitable 
class of moments, a function T : A -> F must be specified such that J 7 realizes the moments in A, 
i.e. ( mJ 7 (p )) = p for all p G A, and F{(mf)) constitutes a suitable (in a sense to be made more 
precise below) approximation to the solution / of the Boltzmann equation (JT]) . Approximating the 
moments in cl by p ~ (m/) and replacing / in (11411 by the approximation J-(p), one obtains the 
following closed system for the approximate moments: 

d t p + d Xi {v i mJ 7 (p)) = (mC(J(p))). (16) 

The closed moment system (1161) is essentially defined by the polynomial subspace, M, and the clo¬ 
sure relation, J 7 . A subspace/closure-relation pair (M, J 7 ) is suitable if the corresponding moment 
system (Ha is well posed and retains the fundamental structural properties of the Boltzmann equa¬ 
tion HD as described in section[2l viz., conservation of mass, momentum and energy, Galilean invari¬ 
ance and dissipation of an entropy functional. Auxiliary conditions may be taken into consideration, 
e.g. that the fluxes and production terms can be efficiently evaluated by numerical quadrature. 

It is noteworthy that moment systems can alternatively be conceived of as Galerkin subspace- 
approximations of the Boltzmann equation in renormalized form. This Galerkin-approximation inter¬ 
pretation can for instance prove useful in constructing error estimates for (|16l) and in deriving struc¬ 
tural properties. In addition, the Galerkin-approximation interpretation conveys that smooth func¬ 
tionals of approximate distributions obtained from moment systems, such as velocity moments, gen¬ 
erally display superconvergence under hierarchical-rank refinement, in accordance with the Babuska- 
Miller theorem; see [2] and also Section HJ Consider the subspace M and let /? : M —y F denote a 
renormalization map. Denoting by V((0, T) x 12; M) a suitable class of functions from (0, T) x 12 
into M, the moment system HHD can be recast into the Galerkin form: 

Find g G V((0,T) x 12; M)) : 

(■ md t /3(g)) + ( mVid Xi /3(g )) = < mC{/3{g ))) Vra G M, a.e. (t,x) G (0 ,T) x 12. (17) 

To elucidate the relation between m and HZD, we associate to f3 : M —» F a function Tp : @{J-p) —> 
F such that J-p(p) — /3(g p ) with g p according to ( mf3(g p )) = p. The domain 9>{J-p) is implicitly 
restricted to moments p G that can be realized by some g G M. The equivalence between the 
Galerkin formulation H3 and the moment system na now follows immediately by noting that 
{rrijI'Aj constitutes a basis of M and inserting g p for g in (1171) . 

In the remainder of this section we review the celebrated moment closures of Levermore (33] and 
Grad [TS] to provide a basis for the subsequent divergence-based moment closures in section [I] 
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3.1 Levermore’s Entropy-Based Moment Closure 

The moment-closure relation of Levermore [33] is essentially characterized by the renormalization 
map /?(•)= exp(-). For this closure relation, a subspace M is considered to be admissible if it satisfies: 

1) ECM; 

2 ) M is invariant under the actions of T u and To] 

3) M c := {m € M : (exp(m)) < oo} has a nonempty interior in M. 

The first condition insists that M contains the collision invariants, which ensures that the moment 
system imposes conservation of mass, momentum and energy. These conservation laws must be 
obeyed if any fluid-dynamical approximation is to be recovered. The second condition dictates that for 
all m £ M, all u £ M. D and all orthogonal tensors O it holds that m{u — (•)) £ M and m(0*(-)) £ M. 
This condition ensures that the moment system exhibits Galilean invariance. As argued by Junk [28] , 
rotation and translation invariant finite dimensional spaces are necessarily composed of multivariate 
polynomials. The third condition requires that M contains functions m such that (3 (m(-)) is Lebesgue 
integrable on R 33 . For /?(•) = exp(-) and M composed of multivariate polynomials, this condition 
implies that the highest-order terms in any variable in M must be of even order. The subset M c 
then corresponds to a convex cone, consisting of all polynomials in M for which the highest-order 
terms in any variable are of even order and have a negative coefficient. One can infer that exp(-) 
maps M c to distributions with bounded moments and fluxes, i.e. g £ M c implies \(ml3(g))\ < oo and 
\vm/3(g))\ < oo for all m £ M. 

In [32] the moment-closure relation associated with /?(■) = exp(-) is derived by minimization 
of the entropy with density i%(f) := /log / — /, subject to the moment constraint. Specifically, 
considering any admissible subspace M, Levermore formally defines the closure relation p i—>• T\fp) 
according to: 

T\.{p) := arg min {(/log / — /) : ( mf) = p} (18) 

/€F 

To elucidate the fundamental properties of the closure relation m, we consider an admissible 
subspace M and we denote by D the collection of all / £ F that yield moments p = (mf) for which 
the minimizer in m exists. The operator J r L ((m(-))) : D —> I is idempotent and its image I C D 
admits a finite-dimensional characterization. In particular, it holds that log I coincides with the 
convex cone M c . The idempotence of the operator Th((m(-))) and its injectivity in D imply that (1181) 
corresponds to a projection. This projection is generally referred to as the entropic projection [22 ] . 
A second characterization of (fT51) follows from the following sequence of identities, which holds for 
any / £ F such that (mf) — p and all Maxwellian distributions, M. — exp(^;) with £ E: 

(/log(//AJ)) = (/(log/-?/’)) = (/log/ - /) + (/(l — , 0)> = (/log/ - /) +« • P (19) 

for some at £ R M . Noting that ot ■ p is independent of /, one can infer from (a and m that Tj, 
according to m is the distribution in F that is closest to equilibrium in the Kullback-Leibler 
divergence, subject to the condition that its moments (m(-)) coincide with p. Similarly, it can be 
shown that T\, according to m minimizes (/ log /) subject to (mf) = p. Therefore, the information 
interpretation of the entropy (/log/) (see Section [2]) enables a third characterization of (1181) . viz. as 
the least-biased distribution given the information (m(-)) — p on the moments. 

The exponential form of the renormalization map associated with csd can be derived straight¬ 
forwardly by means of the Lagrange multiplier method. Provided it exists, the minimizer of the 
constrained minimization problem (1181) corresponds to a stationary point of the Lagrangian (/, a) i —> 
(f log f — f) + at ■ (p — (mf)). The stationarity condition implies that log / a ■ m vanishes, which 
conveys the exponential form / = exp(a • m). It is to be noted that the Lagrange multipliers have 
to comply with an admissibility condition related to integrability. In particular, a ■ m must belong 
to the convex cone M c . 

In [23] it is shown that the moment system (fTOl) with closure T\. corresponds to a quasi-linear 
symmetric hyperbolic system for the Lagrange multipliers. Application of the chain rule to © with 
T\Xp) = exp(a • m) (with, implicitly, p — ( mexp(ct • m))) yields: 

dot \ - , . dot 

A °^~m + z^ Ai ^d^ = s ^ 

i= 1 1 


( 20 ) 



with Ao(a) = (m®mexp(ot ■ m)), Ai(a) — (vim®mexp(a ■ m)) and s(a) = (mC(exp(cc • m))). 
The symmetry of Ai (i — 0,1,..., D) and the positive definiteness of Aq are evident. By virtue of its 
quasi-linear symmetric hyperbolicity, the system (1201) is at least linearly well posed [33]. Moreover, 
under auxiliary conditions on the initial data, local-in-time existence of solutions can be established; 
see, for instance, m- 

Levermore’s moment systems retain the fundamental structural properties of the Boltzmann 
equation. The conservation properties and Galilean invariance are direct consequences of conditions 1. 
and 2. on the admissible subspaces, respectively. Dissipation of the entropy can be inferred 

from the Galerkin formulation ed. by noting that for /}(•) = exp(-) it holds that log/3(-) : M —> M. 
Hence, if g complies with ED and /?(•) = exp(-) then the following identity holds on account of 
Galerkin orthogonality: 

<log (3(g)d t (3(g)) + (log P{g)vid x M) = <log/?( 5 ) C(/?(s))> (21) 

The left-hand side of this identity coincides with d t (r]h(l3(g))) + d Xi (viT]i,(j3(g))), while the right-hand 
side equals (C(/3(g)) rj' L ((3(g))). For g according to (fI71) . the distribution /3(g) — exp (g) thus obeys 
the entropy dissipation relation © with entropy density r] l- 

Levermore’s consideration of entropy-based moment-closure systems in [33], as well as the above 
exposition, implicitly rely on existence of a solution to the moment-constrained entropy minimization 
problem (1181) . It was however shown by Junk in the series of papers [5S][571[55] that for super¬ 
quadratic M the closure relation (fl8l) is impaired by non-realizability, i.e. a minimizer of (1181) may 
be non-existent. Moreover, the class of local equilibrium distributions generally lies on the boundary 
of the set of degenerate densities. In [26], Junk also establishes that the flux ( Vim(3(g)) can become 
unbounded in the vicinity of equilibrium, thus compromising well-posedness of (1201) . The singularity 
of the fluxes moreover represents a severe complication for numerical approximation methods; see 
also [39] . 

The realizability problem of Levermore’s entropy-based moment closure has been extensively 
investigated; see, in particular, [5911571155112311491133] . In [5511551135] it has been shown that the set 
of degenerate densities is empty if and only if the set {a € : mexp(a • m) E L 1 (M' D ,R m )} 

of Lagrange multipliers associated with integrable distributions is open. This result implies that 
degenerate densities are unavoidable for super-quadratic polynomial spaces, because the Lagrange 
multipliers associated with equilibrium are then located on the boundary of the above set; see 
also [23] ■ To bypass the realizability problem, Schneider [39] and Pavan 23 considered the following 
relaxation of the constrained entropy-minimization problem: 

arg min{(/ log / - /) : (mf) <* p} (22) 

/€F 

where the binary relation <* connotes that the highest order moments of the left member are 
bounded by the corresponding moments of the right member. The relaxation of the highest-order- 
mornent constraints serves to accommodate that minimizing sequences {/„} C F subject to the 
constraint (mf n ) — p converge (in the topology of absolutely integrable functions) to an exponential 
density with inferior liighest-order moments; see [56ll5Tll591l23U43j . The analyses in [391143] convey that 
the relaxed minimization problem indeed admits a unique solution, corresponding to an exponential 
distribution. The exponential closure can therefore be retained if the closure relation is defined 
by (1221) instead of (1181) . It is to be noted however that the closure relation (1221) does not generally 
provide a Injection between the Lagrange multipliers and the moments. Moreover, the aforementioned 
singularity of fluxes near equilibrium is also inherent to (1221) . 

Another formidable obstruction to the implementation of numerical approximations of Lever¬ 
more’s moment-closure systems are the exponential integrals that appear in (fT6l) . The evaluation 
of moments of exponentials of super-quadratic polynomials is generally accepted to be intractable, 
and accurate approximation of such moments is algorithmically complicated and computationally 
intensive; see, in particular, [351 Sec. 12.2] and [5B] Sec. 6]. 


3.2 Grad’s Hermite-Based Moment Closure 

In his seminal paper [19] . Grad proposed a moment-closure relation based on a factorization of the 
one-particle marginal in a Maxwellian distribution and a term expanded in Hermite polynomials; see 
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also [201 Sec. V]. The expansion considered by Grad writes: 




■MM EE 

fc =0 i k 


k\ 




(23) 


where c denotes peculiar velocity, if. — (ii,* 2 , ■ ■ ■, *k) is a multi-index with sub-indices fi.) 6 
{1, 2,..., Z)}, a^' are the polynomial expansion coefficients and are D-variate Hermite poly¬ 

nomials of degree k: 


jpW {x ) = (~l) fc d k u(\x\) 

%k u(\x\) dxi t dxi 2 ■ ■ ■ dx ik 


with w(s) 


1 

(2n) d / 2 


exp(— s 2 / 2 ). 


(24) 


The Maxwellian in (1231) can either correspond to a prescribed local or global Maxwellian, or it 
can form part of the approximation; see [191120] . In the latter case, the coefficients associated with 
invariant moments are fixed and it holds that a® = 1, — 0 and a^P — 1. By virtue of the 

specific properties of Hermite polynomials, moments of Grad’s approximate distribution (l23l) can be 
evaluated in closed-form. 

The linear hull of the Hermite polynomials { }o<fc<n coincides with the class of D-variate 
polynomials of degree at most n. The Hermite polynomials in (1241) do not provide a basis of the 
polynomials, however, on account of linear dependence; evidently, the Hermite polynomial in (1241) is 
invariant under permutations of its indices. In [19] . uniqueness of the coefficients in (l23l) is restored 
by imposing auxiliary symmetry conditions on the coefficients. 

Grad’s moment systems can be conveniently conceived of as Galerkin approximations of the 
Boltzmann equation in renormalized form in accordance with (ED. For a prescribed Maxwellian, the 
renormalization map simply corresponds to /3 : < 71 —> A4g. Incorporation of the Maxwellian in (1231) in 
the approximation can be represented by the renormalization map: 


/3 : 3 exp(77 E g) x (l + (Id - II E )g) 


(25) 


where 7Z E :Mh)E denotes the orthogonal projection onto the space of collision invariants and Id 
represents the identity operator. The embedding ECM implies that 77 E M = E and (Id — 77 E )M = 
M\E. Hence, the projection in (1251) provides a separation of M into E and its orthogonal complement. 
It is notable that the renormalization map in Grad’s moment system can be conceived of as a 
linearization of Levermore’s exponential closure relation in the vicinity of A4. In particular, setting 
ip — logTW € E, the following identities hold pointwise: 

exp(fif) = exp(ip) exp (5 - ip) = Mex.p{g-ip) = M (l + (fir - ip) + 0(\g - ip\ 2 )) (26) 

as (g — ip) —>• 0. To derive the renormalization map j3 : g 1 —>■ M.g for prescribed Maxwellians, it 
suffices to note that 1 + g — ip € M. To infer the renormalization map (1251) if A4 is retained in the 
approximation, we note that setting ip — II E g and omitting the remainder in (l26|) yields (1251) . 

For a prescribed (global or local) Maxwellian A4, Grad’s moment systems dissipate the entropy 
rj x ^{f) ■— ^A4(//A4 — l) 2 , provided that r ] x 2 represents an entropy density for the collision operator 
under consideration. It can for example be shown that r ] x 2 is generally a suitable entropy density 
for collision operators linearized about A4 (see ED) and for BGK collision operators. Dissipation 
of the entropy ( 77 ^ 2 ) can be directly inferred from the Galerkin formulation (fT71) . by noting that for 
f3(g) — Mg it holds that: 

r! x ., (£(<?)) = 0{g)/M - 1 = g - 1 G M (27) 

Hence, rp 2 resides in the test space M in ED and dissipation of (rj x 2 ) follows from Galerkin orthog¬ 
onality. 'f'he entropy ( 77^2 (/)) can be associated with the ip x 2 -divergence of / relative to A4 with 
<p x a(s) — ^(s — l) 2 ; cf. (HID . Grad’s moment-closure relation can in fact be obtained by minimization 
of the <p x 2 -divergence subject to the moment constraints: 


F G (p) = arg min {(//A4)) : (mf) = p} (28) 

fe f 


The minimization problem (1281) is not impaired by the realizability problem inherent to H 8 D , because 
the moment functionals (m(-)) are continuous in the topology corresponding to (g x 2 ). 
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If the Maxwellian is retained in the approximation, then an entropy for the corresponding moment 
systems can be non-existent or its derivation is intractable. However, for any entropy density g for 
the collision operator, the following identity holds by virtue of the Galerkin-orthogonality property 
of /3 := /3(g) in (HU): 

dt(d(P)) + d Xi {vMP)) = ((v'(P) ~ m)d t p) + ((g'(P) - m)vid Xi (3) 

for arbitrary me M. Equation (l29l) implies that solutions to Grad’s moment systems dissipate any 
entropy (g) for the collision operator up to inf mG M || g'(/3(g)) — m||, in some suitable norm || • ||. For 
example, introducing the condensed notation go — Il^g, gi — (Id — I7e)<? and the convex functional 
?? : E x M \ E —» R. according to 

v(9o,9i) = {go - l)e So (l + gi) + e 9o gi(l + gi) (30) 

the renormalization map in (1251) corresponds to /3(g) — e So (l + 51 ) and it holds that 

dv(go,gi) = (go + gi)e 90 (1 + gi) dg 0 + (g 0 + 2gi)e 90 dgi 

= (go + gi)(d go /3dg 0 + d gi /3 dg x ) + g ie so dg 1 (31) 

= (go + gi) d/3 + gie 90 dg 1 

Considering that go + gi e M, it follows from (1291) that if g in (1301) is an entropy density for the 
collision operator, then Grad’s moment systems with /?(•) according to (l25l) dissipates g up to 0(g\) 
as g\ vanishes (in some appropriate norm). Note that g\ vanishes at equilibrium. 

Grad’s moment-closure relation exhibits several fundamental deficiencies that may cause break¬ 
down of the physical and mathematical structure of the corresponding moment-closure system for 
large deviations from equilibrium. First, the expansion m admits inadmissible, locally negative 
distributions. Second, the moment systems are generally non-symmetric and hyperbolicity is not 
guaranteed. It has been observed in mm that Grad’s moment-closure systems can indeed exhibit 
complex characteristics and loss of hyperbolicity. 


4 Divergence-Based Moment Closures 

In this section we present a novel moment-closure relation based on an approximation of the ex¬ 
ponential function. The considered approximation is derived from truncations of the standard limit 
definition of the exponential exp(-) := lim„_ >OQ (l + (■)/n) n ss (1 + (-)/N) N . It is noteworthy that 
unlike the exponential function, in the limit as v —> —00 the truncated exponential as well as its 
derivative do not vanish. The former condition is needed to preserve the decay properties of the 
exponential function while the latter condition is needed to preserve the same absolute maximum 
and minimum as the exponential. Moreover, as opposed to the exponential function, the truncated 
exponential can be negative if N is odd. Several approximations of the exponential function that 
preserve the aforementioned properties of the exponential have been proposed in the literature; see, 
for example, [5511421129] and references therein. These so-called deformed exponentials can generally 
serve to construct moment-closure renormalization maps, with properties depending on the particular 
form of the deformed exponential and the construction. In [55], Tsallis proposed the q- exponential: 

exp q (x) := (l + (1 - q)x)+ (1 ~ q) (32) 

with q 7 ^ 1 and (•)+ = |(*) + \ H the non-negative part of a function extended by 0. The g-exponential 
in (E3) is related to the non-negative part of the truncated limit definition of the exponential by 
1 — q — 1 /N . We will consider renormalization maps of the form 

Pn ■ 9 ^ Mexp q (g) = M(l + -^) + (33) 

with A4 a prescribed local or global Maxwell distribution. The renormalization map /3n can be 
construed as an approximation to the exponential renormalization map about the Maxwellian dis¬ 
tribution A4. We will establish that the moment-closure distribution m can be derived as the 
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minimizer of a modified entropy that approximates the Kullback- Leibler divergence near M and 
that belongs to the class of ip- divergences. In addition, we will show that the resulting moment sys¬ 
tem overcomes the aforementioned deficiencies of Grad’s and Levermore’s moment systems, while 
retaining the fundamental properties of the Boltzmann equation presented in Section [21 
The renormalization map (13311 engenders the following moment-closure relation: 

( r\ . YY1 \ N 

1 + n ) ( 34 ) 

where the moment densities p and the coefficients a. are related by p = (m A4 exp (a ■ m)). Given 
a polynomial subspace M D E with a Galilean-group property (admissibility conditions 1 and 2 
in section ELD, the moment system corresponding to fTTTTIi conforms to (fl6l) with, in particular, the 
moment-closure relation T N according to (1551) . 

To elucidate some of the characteristics of the renormalization map (1331) , we regard it in compari¬ 
son with the renormalization maps associated with Levermore’s exponential moment-closure relation 
and Grad’s moment-closure relation with a prescribed Maxwellian prefactor. The renormalization 
map associated with Levermore’s moment-closure relation is given by g exp(g); see Section l3Tl By 
virtue of the vector-space structure of M D E, for an arbitrary Maxwellian distribution A4 it holds 
that log Af + M = M. Hence, for g E M, the renormalization map g H > exp (g) can be equivalently 
expressed a,s g ^ exp(log Af + g). In the limit N —> oo, we obtain for (1331) : 

lim fi N (g) — M lim + = M expfo) = exp(log M + g) (35) 

TV—»oo TV—,oo V A / + 

Equation (1351) implies that in the limit N —> oo, the renormalization map in (1331) coincides with the 
exponential renormalization map associated with Levermore’s moment-closure relation. For finite N, 
the moments ( mPi\r(g )) and fluxes (mvf}N{g)) with m,g E M correspond to piecewise-polynomial 
moments of the Gaussian distribution Ad. The evaluation of such moments is tractable, as opposed to 
the evaluation of moments and fluxes for the exponential renormalization map. In addition, for super¬ 
quadratic approximations M D E, the exponential renormalization map associated with Levermore’s 
closure can lead to singular moments and fluxes in the vicinity of equilibrium, i.e. as g approaches 
E. The fundamental underlying problem is the realizability problem; see Section [3.11 and [26] . Ac¬ 
cordingly, one can form sequences {<?„} such that exp(< 7 „) —>• E (in the L 1 topology) while there 
exist m E M such that |(mexp(< 7 „))| —»• oo or |(mwexp(< 7 ri ))| —>• oo. One can infer that due to the 
exponential decay of the prefactor Af and the polynomial form of the renormalization map in (1331) , 
moments and fluxes corresponding to (1551) are non-singular near equilibrium. To compare (1331) to the 
renormalization map corresponding to Grad’s moment-closure relation with a prescribed Maxwellian 
prefactor, g i—>- A 4g (see Section [3721) . we note that by virtue of the vector-space structure of M D E, 
it holds that 1 + M = M. Hence, for g E M, the renormalization map g A ig can be equivalently 
expressed as /3q : g >->• Ad(l + g). Comparison of and /3q imparts that /3± — (/3q) + , i.e. for N — 1 
the renormalization map /?n in (1551) coincides with the non-negative part of the renormalization 
in Grad’s closure, extended by zero. Therefore, the renormalization map (1331) avoids the potential 
negativity of the approximate distribution inherent to Grad’s closure and the corresponding loss of 
hyperbolicity of the moment system. 

The moment system corresponding to (1331) retains conservation of mass, momentum and energy 
as well as Galilean invariance. The conservation properties can be directly deduced from the Galerkin 
form m of the moment system, by noting that E is contained in the test space M, in accordance 
with admissibility condition 1 in section ELD- Galilean invariance is an immediate consequence of 
admissibility condition 2. However, contrary to Levermore’s moment system, the moment system with 
renormalization map (1551) does not generally dissipate the relative entropy (/log (//Ad)), because 
the inverse of /3jy(-) does not correspond to log(-) and, therefore, log/3jv(ff) does not generally belong 
to the test space M for g E M; cf. Section 13.11 The moment system closed by (1551) does however 
dissipate a modified entropy. To determine a suitable entropy function for the moment system with 
renormalization map (1331) . we observe that: 

= — N =: logjy((•)/A / l) 


(36) 
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provides an inverse of /3 jv according to (l33l) with domain M>o- The function log yields an approx¬ 
imation to the natural logarithm, corresponding to the inverse of the q-exponential in (1321) . The 
approximate logarithm is eligible as the derivative of an entropy density associated with the moment 
system with renormalization (1331) . In particular, defining the entropy density as 

/ N 2 ( f \ 1/N \ N 

= ~ n ) +m ttn < 37 > 

it holds that r/' N — /dj^ 1 2 3 and, hence, v'n(Pn(-)) ■ M —> M. The constant in (l37l) has been selected 
such that 77 at(AI) vanishes. The entropy corresponding to (1371) can be cast in the form of a relative 
entropy associated with a (^-divergence, in accordance with m ■ To this end, we introduce 

™(') = (')(rSv ( ' )1/w - JV ) + Tfjv < 38 > 

and note that J7 at(/) = M.ip N (f / M). Convexity of the function (ppf and of the corresponding entropy 
density T)l follows by direct computation: 

M) = (') 1/JV - 1 (39) 

Therefore, ip" N is strictly positive on R>o- Moreover, it holds that <£>jv(l) = 0. In conclusion, if 
tjn{-) — M<p((-)/M) is an entropy density for the collision operator C according to (f5j), then the 
approximate distribution /3N(g) of the moment system (fT71) with renormalization map (l33l) complies 
with the local entropy-dissipation relation: 

d t (vN{PN(g ))} + d Xi {vir] N (p N (g))) = (C(/3 N (g))ri' N (p N (g))) < 0 (40) 

We recall that the premise on the collision operator is for example satisfied by the BGK and extended 
BGK operators; see also Remark [1] 

The moment-closure relation (1341) can be derived by minimization of the (^Ar-divergence subject to 
the moment constraint; cf. the definition of Levermore’s closure relation according to (1181) . Consider 
the constrained minimization problem: 

T N (p) := arg min { (t7at(/)) : (mf) = p} (41) 

fe f 

Formally, the solution to m can be obtained by the method of Lagrange multipliers. The minimizer 
in (1411) corresponds to a stationary point of the Lagrangian (/, a ) (r/]\r(f)) + a ■ (p— (mf)). The 

stationarity condition implies that r/' N (f) — a-m — 0 and, on account of (1361) . that /3 a) 1 (/) = am. It 
follows directly that the minimizer in (1411) is of the form ^^(p) — j3^(a ■ m) in conformity with (1341) . 
Contrary to the entropy minimization problem (1181) underlying Levermore’s closure relation, the 
minimization problem m is well posed. Existence of a solution to the minimization problem m 
can be deduced from results for generalized projections for non-negative functions by Csiszar in CH¬ 
IU mi it is shown that the minimization problem 

inf | fMvMm/m) v(dv) : /i e x| (42) 

over a constrained set of non-negative functions, 

X= j/: J aj(v)f(v)v(dv) = pj, j € /j| (43) 

for certain countable functions possesses a minimizer belonging to X provided that the 

following (sufficient) conditions hold: 

1) X is a convex set of non-negative functions and the infimum in (l42l) is finite; 

2) lim^oo tp'(s) = oo; 

3) J t P*(£\aj(v)\)f 2 (v) v(dv) is finite for all £ > 0 and j E IJ. 
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The function ip* in condition [3] corresponds to the convex conjugate of p. Comparison conveys 
that (TO1) conforms to (l42l) - (l43l) with fa — M, z'(-) Lebesgue measure and {aj}j 6 /j a monomial 
basis of M. Convexity of the constrained distributions follows from the linearity of the moment 
constraints. Finiteness of the infimum is ensured by the fact that the inhmurn over the constrained set 
is bounded from below by the infimum over the unconstrained set, and the latter attains its minimum 
of 0 for f — JA. The minimization problem (Tfl1) thus complies with condition |TJ Compliance with 
condition [2] follows from p' N (s) — Ns 1 ' N — N and lim^oo s 1//Ar = oo. To verify condition[3j we note 
that the convex conjugate of ipu is: 

N (( t \ N+1 \ 

v " {t) = ^J st - VN{s)) = TTN\{ 1 + N) + - 1 ) <44) 

Condition [5] therefore translates into the requirement that 

(M V S f (fK|)> = T ^C^(l + ^l)" +1 -A(^ (45) 

is bounded. By virtue of the exponential decay of the prefactor A4 and the fact that | m,j \ N+1 
increases only polynomially, the expressions in (l45l) are indeed finite for any £ > 0. The minimization 
problem (Tfl1) therefore also satisfies condition [3] It is notable that the minimization problem (1181) 
pertaining to Levermore’s moment closure satisfies conditions [U and [5] but violates condition Q5J 
To establish that the closure relation (l34l) leads to a symmetric-hyperbolic system, we insert El 
into the generic form m of moment systems, and note that application of the chain rule and product 
rule leads to a system of the form (12011 with: 


A 0 {a) = (m®mM (l + + ^ 

/ / . yyi \ N — l \ 

Ai(a) — (vim ®mM -J + ) 


, , , . . „ / a ■ m\ N 

s(a) = (mClM (l + 


-(ml 


N 


fdM dM\ 

-) + Ur + 2 >&rJ 

i=i * 


The symmetry of Aq, ..., Ap is evident. Positive-dehniteness of Aq follows from: 


, ot-m\ N - 1 \ , . 9 / am\ N ~ . 

7 '(m« mM (^1 H-J ) 7 = ((7 • m) 2 M ^1 H-J )>0 


(46a) 

(46b) 

(46c) 


(47) 


The inequality in (1471) reduces to an equality if and only if 7 = 0 or am — —N. The latter case is 
pathological, because am — —N implies that (p) = 0. 

It is noteworthy that the second constituent of the production term s(a), i.e. the term repre¬ 
senting the contribution of d t M. + Vid Xi M to the production, may cause blow up of solutions to 
the hyperbolic system m with (l46l) in the limit t —> 00 . Hence, the hyperbolic character of (1161) 
with (1461) ensures stability of solutions only in finite time. If M corresponds to a global Maxwellian, 
then d t A4+Vid Xi A4 vanishes and the stability provided by hyperbolicity also holds in the ad-infinitum 
limit. 


5 Numerical results for the ID spatially homogeneous Boltzmann-BGK equation 

To illustrate the properties of the moment-system approximation (El) with the divergence-based clo¬ 
sure relation encoded by the renormalization map (1331) . this section presents numerical computations 
for the spatially homogeneous Boltzmann-BGK equation in ID: 

dtf = -r- 1 (f-£ f ) 

/(0,u) = fo(v) 


(48a) 

(48b) 
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for some given initial distribution /o- The corresponding moment system writes: 

(49a) 
(49b) 


d t (mJ 7 N ) = -r 1 <m(J r /v - £t n )) 
Fn(0, v) = (.Fjv)o(w) 


with Tjy according to (1M1) and (J 7 at)o defined by the minimization problem (TUI) subject to the 
moments corresponding to the initial distribution: 

(Fn )o :=argmin{(r?Ar(/)) : {mf) = (mf 0 )} 

f GV (50) 

={/ E¥ : f — /3 N (g),g G M, {mf) = (m/ 0 )} 


The systems in (1481) and m represent initial-value problems for the ordinary differential equa¬ 
tions (148al) and (149al) . The solutions of the initial-value problems (1481) and (1491) are, respectively, 


f{t, v) = e - ( t-t °)/ r / 0 + (1 - (51a) 

= e-^/ T (T N ) 0 + (1 - e-^/ T )£ TN (51b) 

For the considered spatially-homogeneous case, the collision-invariance properties of the collision 
operator imply that £f — £f 0 and, similarly, £t n — £(j r N ) 0 - Furthermore, the constraints in the 
minimization problem in (1501) impose £f 0 — £(f n ) 0 - Based on the expressions for the solutions in (1511) . 
it then follows that: 

||/(M- -MV)|| =e-( t - to)/r ||/ 0 -(^)o|| (52) 

in any suitable norm. Equation (!52l) conveys that the accuracy of the approximation of f(t,-) by 
•) at any time t > 0 depends exclusively on the accuracy of the approximation of the initial 
condition /o by (J-n)o according to (1501) . In the remainder of this section we therefore restrict our 
considerations to numerical examples that illustrate the approximation properties of ^w-divergence 
minimizers and to properties of the projection problem (1501) . 

We consider approximations of the distributions: 

e -h(y-2f e ~h( v + 2 ) 2 

fliv) = + sftn ’ 

e -U v ~ 2 f e-J^+ 2) 2 
h{v) = + ’ 

-2(v-2) 2 -§v 2 e ~(v+2) 2 

h{v) = -j7jT + 7 ^ 74 + 

by means of moment-constrained <^]\r-divergence minimizers in polynomial spaces of increasing or¬ 
der. The distributions in (l53aD - (l53cjl . shown in Figure [TJ correspond to distributions of increasing 
complexity, viz., a symmetric bi-modal distribution, a non-symmetric bi-modal distribution and a 
non-symmetric tri-modal distribution, respectively. For the pre-factor A4 in the renormalization 
map (l33l) and, accordingly, in the relative entropy (t]N(f)) = {M.g>N{f / M)) associated with the 
(/jjv-divergence in (l38l) . we select the global equilibrium distribution £f 0 . I 11 particular, denoting by 
M/c = span{l, v ,..., u fc-1 } the space of polynomials of degree k— 1, the minimization problem in (1501) 
with k moment {k > 3) constraints engenders the nonlinear-projection problem: 

(Jqv)o = M (l + jj) + , g G M fc : = ( m fo) Vm e M k (54) 

Expanding g{v) — onv l ~ x — a. ■ m(v), Equation (1541) corresponds to a nonlinear algebraic system for 
the coefficients a. To evaluate the integrals in (l49l) . we first determine the roots of the polynomial 
(1 + a. ■ m/N ) and then establish the limits of the positive parts to compute the corresponding 


(53a) 

(53b) 

(53c) 
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Fig. 1: Distributions fi (left), /2 ( center ) and (right) according to (1531) and the corresponding 
approximations with k — 3 (dotted), k — 6 (dash-dot) and k — 15 (dashed) moments obtained from 
the moment-constrained (^-divergence minimization problem (1501) . 


contributions to the moments of J\A(1 + a ■ m/N)f_. The integrals are evaluated by applying a 
suitable transformation of the integration variable and invoking the following rule: 


v k dv = 


1 


r 


~2 lr 


1 + k 
2 

1 + k 


-r 


- r 


1 + k 
2 ,V 
1 + k 


sign 1 +fc (ui) 


(55) 


-y 0 ) ) si gn 1 +/£ (u 0 ) 


1 +k/ 


where —oo < vq < oo and —oo < v\ < oo and r(-) and F(-,-) are the complete and incomplete 
gamma functions respectively. The coefficients a are extracted from the system (El by means of 
the Newton method. It is to be noted that (1 + g/N)^_ is Frechet differentiable with respect to g by 
virtue of the fact that, evidently, changes in the sign of 1 + g/N occur only at roots. A consistent 
Jacobian for the tangent problems in the Newton method is provided by: 



The Jacobian matrix J(a) in the right member of (1551) can be identified as a symmetric-positive 
definite matrix and, hence, the tangent problems in the Newton method are well posed. The Jacobian 
is however of Hankel-type and it becomes increasingly ill-conditioned as the number of moments 
increases; see for example mm- Consequently, the convergence behavior of the Newton process 
deteriorates for higher-moment systems. To illustrate the dependence of the convergence behavior 
of the Newton process on the number of moments, Figure [2] (left) plots the ratio of the 2-norm 
of the update in the Newton process, ||<5 ck|| 2 , over the 2-norm of the solution vector, |||| 2 , 
versus the number of iterations for polynomial orders k — 7, 9, 11, 13 for the three test distributions 
in m • The ratio ||Ack|| 2 /1|ck^ 77 ' -1-1 ^ ||2 can be conceived of as the relative magnitude of the update 
vector. Figure[2] (right) plots an approximation of the corresponding infinity-norm condition numbers, 
x 00 (o:) = ||J(a)|| 00 || J _ 1 (o:)|| 00 , of the Jacobian matrices. The results in Figure [2] convey that 
the condition number increases significantly as the number of moments increases. For k — 7 the 
condition number is approximately 10 3 , while for k — 13 the condition number exceeds 10 5 and can 
even reach 10 10 . For high-order approximations, the convergence behavior of the Newton process is 
generally slow and non-monotonous. However, in all cases the relative update can be reduced to a 
tolerance of 10 _4 . 

To illustrate the approximation properties of the moment method with closure relation (1341) , Fig¬ 
ure [3] (left) plots the L l (M)-norm of the relative error in the approximation 7F\ to the distribution fi 
(i — 1,2,3) according to (l53l) and Figure [3] (right) the corresponding relative error in the cosine 
moment, 


Wfi . = |(cQs(-)/i) - (cps(-)j r 2 )| 

Wfi 11L 1 (M) ’ 0112 l(cos(-)/;}| 


(57) 


respectively. The cosine moment serves to investigate the super-convergence properties of the ap¬ 
proximation in accordance with the Babuska-Miller theorem |2] ; see Section [3] A non-polynomial 
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Fig. 2: Convergence of the Newton process for the nonlinear projection problem (15411 and condi¬ 
tioning of the corresponding Jacobian matrices: (left) relative magnitude of the Newton update, 
|| Ja|| 2 /||a:( rl+1 ) || 2 , versus the number of iterations for k — 7, 9,11,13 and for distributions (153aD (top), 
(I53bl) (center) and (I53cl) (bottom)-, (right) corresponding oo-norm condition numbers, >c 00 (a^ n '), of 
the Jacobian matrices according to (1561) . 


moment has been selected to examine the convergence behavior, because for any polynomial moment 
(off) with <t E M; the approximation (oJ-lf) provided by the fc-moment approximation is exact 
for all k > l, on account of the constraints in (1501) . Figure [0] (left) indicates that ||/j — J^IUmk) 
converges exponentially with increasing k, i.e., there exist positive constants C and £ such that 
||fi — < C(f~ k . In particular, £ ps 10 -013 0.74 for the bi-modal distributions f\ and fa 

and C ps io -0 ' 085 ps 0.82 for the tri-modal distribution fa. Comparison of the left and right panels 
in Figure [3] conveys that the approximation of the cosine moment indeed converges at a higher rate 
than the L 1 (M)-norm of the approximation itself. Figure [3] (right) conveys that the cosine moment 
converges at a rate of £ ps IO -0 " 58 ps 0.26 for both the bi-modal distributions fa, fa and the tri-modal 
distribution fa. 


6 Conclusion 

To avoid the realizability problem inherent to the maximum-entropy closure relation for moment- 
system approximations of the Boltzmann equation, we proposed a class of new closure relations 
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Fig. 3: Approximation properties of the moment method with closure relation (1341) : (left) convergence 
of the relative error in the L 1 (M)-norm erri according to (1571) for /i ,/2 and fa in (1531) : (right) 
corresponding convergence of the relative error in the cosine moment, err 2 - 


based on ^-divergence minimization. We established that (^-divergences provide a natural gener¬ 
alization of the usual relative-entropy setting of the moment-closure problem. It was shown that 
minimization of certain (^-divergences leads to suitable closure relations and that the corresponding 
moment-constrained (^-divergence minimization problems are not impaired by the realizability prob¬ 
lem inherent to relative-entropy minimization. Moreover, if the collision operator under consideration 
dissipates a (^-divergence, then the corresponding minimal-divergence moment-closure systems re¬ 
tain the fundamental structural properties of the Boltzmann equation, namely, conservation of mass, 
momentum and energy, Galilean invariance, and dissipation of an entropy, sc. the (^-divergence. For 
suitable (^-divergences, the closure relation yields non-negative approximations of the one-particle 
marginal. Divergence-based moment systems are generally symmetric hyperbolic, which implies lin¬ 
ear well-posedness. 

We inferred that moment systems can alternatively be conceived of as Galerkin approximations of 
a renormalized Boltzmann equation. We considered moment systems based on a renormalization map 
composed of Tsallis’ ^-exponential. This renormalization map is concomitant with a (^-divergence 
corresponding to the anti-derivative of the inverse ^-exponential, which yields a natural approx¬ 
imation to relative entropy. The evaluation of moments of ^-exponential, elementary in numerical 
methods for the corresponding moment system, is tractable, as opposed to the evaluation of moments 
of exponentials of arbitrary-order polynomials, connected with maximum-entropy closure. 

Numerical results have been presented for the one-dimensional spatially homogeneous Boltzmann- 
BGK equation. The nonlinear projection problem associated with the moment-constrained (^-divergence 
minimization problems was solved by means of Newton’s method. We observed that the condition 
number of the Jacobian matrices in the tangent problems generally deteriorates as the number of 
moments increases. Nevertheless, in all considered cases approximations up to at least 14 moments 
could be computed. We observed that the ^-exponential approximation converges exponentially in 
the L (M)-norm with increasing number of moments. Moreover, we demonstrated that functionals 
of the approximate distribution display super convergence, in accordance with the Babuska-Miller 
theorem for Galerkin approximations. 


A Generalized BGK collision operator 

In [33j , Levermore introduced a class of multiscale generalizations of the BGK collision operator based on a finite 
sequence of increasingly constrained entropic projections of the form (1181) . In particular, given an admissible 
space of polynomials M, consider a sequence of nested subspaces with Mo = E and strictly contained 

in M, i.e. E = Mi C M2 C ... C M^ C M. For each k and / £ F, let / F k (f) = : T k be the M^-moments 
constrained entropic projection of /, 

T k (f) := arg min {{ 77 (g)) : (mg) = (mf), Vm £ M*,}, 

g£F 


( 58 ) 
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with g(g) = glogg — g, under the assumption that (1581) admits a solution for each k. Based on the sequence of 
projections {J rfc }^_ 1 , one can define a multiscale relaxation operator: 

K-l 

C(f) = -d K (f -F K )-J2 9 k(F k+1 ~ K k ) (59) 

k =1 

with {9k}^ = i an increasing sequence of positive relaxation rates depending on /. The relaxation rate Ok with 
k 6 {1, 2,. .., K — 1} constitutes the rate at which _F fc + 1 decays to T k , while 6k is the rate at which / decays 
to T K . In |33| it is shown that the Prandtl number can be controlled via the relaxation rates. 

The above construction of the generalized BGK operator can be extended to (/^-divergences. To this end, 
consider an arbitrary (/^-divergence and let / K > F k (f) = '■ T k denote the corresponding divergence-minimization 
projection according to (1581) . i.e. T k is defined by (15811 with g(-) = A4tp((-)/M)■ Based on the projections T k , 
an extended BGK operator can be defined analogous to (1591) . To establish that g'(-) = corresponds 

to an entropy density for the generalized BGK operator, we first note that the (strong) convexity of g implies: 

(g\s)-g\t))(s-t))> 0 (60) 

for all s, t in the domain of 77 and equality in (|60f) holds if and only if s = t. Rearranging the sum in (|59|) yields: 

K — l 

C(f) = -6i (/ - F 1 ) - Y, ( 0fc +i - ° k ) (/ - Tk+1 ) ( 61 ) 

fc=l 

From the minimization problem (15811 we infer that for all k it holds that g'(F k ) £ Mfc and {m(f — F k )) = 0 
for all m £ M*,. Hence, ( g'(F k )f) — (g'(F k )F k ) = 0 yields a partition of zero for all k. From (1611) and the 
aforementioned partition of zero, we obtain 

(v'(f)C(f)) = - t 1 )) - J2 {6k +1 - e h ){n'(j)(f ~ ^ fc+1 )) 

^ K -1 (62) 

= -ei({vV) ~ v'ir 1 )(f - r 1 )) - E ( 9k +i ~ ~ v'(F k+1 )) (/ - -^ fc+1 )) 

fc=i 

From 0i > 0 and 9k +i > 6k (k = 1,. . . , K — l), and the convexity of g(-) according to (1601) we conclude that g and 
C satisfy the dissipation relation (O, i.e. (g'(f)C(f)) < 0 for all admissible /. To verify the second prerequisite 
relation between g and C, viz., the equivalence of the statements in ©, we first observe that the implication 
©(i) => ©(«) is trivial. To validate the reverse implication in ©, we note that ©(ii) in combination with 
the convexity of g according to ([601) and the ultimate expression in (t62l) implies that ( g'(f ) — r?'(-F fc ))(/ — J rk ) 
vanishes almost everywhere for all k = 1,..., K. This, in turn, implies that / = T 1 = • • • = T K . Condition ©(i) 
then follows directly from (ED) . To verify the implication ©(«) => ©) (,.;,), we note that T k according to (1581) 
satisfies g’)T k ) £ M& for all k. Recalling that ©(«) implies / = T 1 , we infer g'(f) £ Mi = E in accordance 
with many Finally, the reverse implication ©(iii) => © (ii) follows immediately from m and the moment 
constraints in (1581) . 
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